%solving strategy using residual based incremental
%update static scheme with combination with
%residual based newton raphson strategy
%in this scheme, lambda is provided in array. This
%helps to perform a prescribed analysis.
%Add in the processInfo which is relevant for
%implicit/explicit integration scheme
function [elements,output] = solving_strategy_static_3...
    (elements,model,lambda,maxits,...
    tol,postProcessInfo)
% tol = 1e-6;
% maxits = 50;
dim = 2;
numTimeStep = size(lambda(:),1);
output = [];
processInfo.integrationType = 'implex';
disp(['perform implicit explicit integration scheme']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time steps iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.preAnalysis();
mesh = model.getMesh();
nDofs = mesh.countDofs();
r = zeros(nDofs,1);
r = model.assembleR(r);
u = zeros(nDofs,1);
for i=1:numTimeStep
    lr = lambda(i)*r;
    if i > 1
        processInfo.Dt = lambda(i) - lambda(i-1);
    else
        processInfo.Dt = lambda(i);
        processInfo.Dt_old = processInfo.Dt;
    end
    
    %update model's states
    model.setX(lambda(i),u);
    [elements,processInfo] = InitializeSolutionStep(elements,model,...
        processInfo);
    
    %assemble ri
    ri = assembleRi(elements,model);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Newton-Raphson iteration
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cnt = 0;
    while 1
        cnt = cnt+1;
        disp(['substep ' num2str(cnt) '--------------']);
        %assemble kt
        kt = assembleKt(elements,model);
        %         rcond(kt)
        %         find(eig(kt) < 0)
        
        du = kt\(lr-ri);
        u = u + du;
        model.setX(lambda(i),u);
        
        %update model's states
        [elements,processInfo] = ...
            CalculateMaterialResponse(...
            elements,model,processInfo);
        
        %assemble ri
        ri = assembleRi(elements,model);
        
        %compute convergent criteria
        %res = norm(du)/norm(u);
        res = norm(r-ri);
        
        disp(['substep ' num2str(cnt) ' complete']);
        if res < tol
            disp(['Newton Raphson iteration converged in ' ...
                num2str(cnt) ' loops']);
            break;
        end
        if cnt>maxits
            disp(['Newton Raphson iteration did not converge in ' ...
                num2str(maxits) ' loops']);
            break;
        end
    end
    
    [elements,processInfo] = FinalizeSolutionStep...
        (elements,model,processInfo);
    
    output = processOutput(i,dim,elements,model,...
        output,postProcessInfo,processInfo);
    
    disp(['time step ' num2str(i) ...
        ' finished with lambda = ' num2str(lambda(i))]);
end

end
